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We present a general formalism for the dissipative dynamics of an arbitrary quantum system in 
the presence of a classical stochastic process. It is applicable to a wide range of physical situations, 
and in particular it can be used for qubit arrays in the presence of classical two-level systems 
(TLS). In this formalism, all decoherence rates appear as eigenvalues of an evolution matrix. Thus 
the method is linear, and the close analogy to Hamiltonian systems opens up a toolbox of well- 
, developed methods such as perturbation theory and mean-field theory. We apply the method to the 

• problem of a single qubit in the presence of TLS that give rise to pure dephasing 1/f noise and solve 

' this problem exactly. The exact solution gives an experimentally observable improvement over the 

^Nj , popular Gaussian approximation. 
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I. INTRODUCTION 



With researchers motivated by the prospect of quantum computing, qubit dephasing has been a topic of intense 
research over the past decade. Various models of this venerable IjJ phenomenon have been investigated. The most 
^ _ popular have been the spin-boson or spin bath models 0, 0, Another important version has been that of an 
• electron spin coupling to nuclei [1, In recent studies of superconducting qubits, however, it has been found that 
^ ! 1 /f-type noise is the chief source of dephasing d, H, [l^j [Hlj [l3 • The sources of this noise are two-level systems (TLS) 
O^. [13I O, [l5| with a wide spectrum of switching rates. It is likely to be important in virtually any solid-state system 
that serves as a host for qubits, as TLS arc ubiquitous in bulk materials. This noise is usually modeled as classical 
noise. 

^ ' Our aim in this paper will be to present a formalism that solves for the dissipative dynamics of an arbitrary quantum 
CO ' system in the presence of a classical stochastic process. This is a very general model of classical noise. The formalism 
\l depends on a combination of the " Liouvillian" app roach to the evolution of the density matrix 16] with methods 
from the classical theory of stochastic processes (l7| |. In particular, the formalism applies to an ensemble of TLS with 
. . any distribution of switching rates and couplings to the quantum system. It has the great advantage of reducing to 
' a linear system of equations, and in fact there is a close analogy to the usual Hamiltonian formulation of quantum 
, mechanics. It is exact, making no approximation as to the strength of the coupling relative to the inverse of the time 
scales of the noise. 

This method has been derived for a specific example in previous work [l8| . As an illustrative case, we use the new 
method to solve the problem of a single qubit in the presence of TLS that give rise to pure dephasing 1/f noise. Other 
solutions of this problem have been found by previous authors [H, [l^] , there have been numerical studies [2II] , and 



the subject has recently been comprehensively reviewed [22[, so this problem is a good testbed for our method. It 
■ also allows us to exhibit the Hamiltonian analogy, which in this case is to a spin 1/2 system. The illustrative case 
points the way to other interesting models that are not exactly solvable, but to which the method also applies. 

This paper is concerned with mathematical methods. Application to specific physical systems will be given in future 
work. The particular case of superconducting qubits has recently been treated [23l| . The main new results of a general 
nature are found in Eqs. 1101 1161 and the physical interpretation following Eq. 1171 New results for strong-coupling 
(1/f and similar) noise are found in Eqs. [Ml and The most convenient starting point for future calculations of the 
effects of 1 /f and other broad-spectrum noise is found in Eq. [50l 
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II. GENERAL METHOD 



We consider the general problem of a quantum system in the presence of classical noise. The quantum system is an 
iV„-state system, so its Hilbert space is 7V„-dimensional. The classical system has Nf. states labeled by the index a. 



The initial state of the composite system is given by the Hermitian iV, 
probability distribution Pa{t = 0). < Pa < 1. p and P satisfy 



q X Nq density matrix p {t 



Tr p{t) = 1 and ^ P^ (t) = 1 



0) and the classical 



(1) 



at all times t. The classical environment passes through a sequence S of discrete states during the course of the time 
evolution. The probability distribution of these states evolves according to the master equation 



dPg jt) 
dt 



J2 ^-bPb (t) . 



(2) 



b=l 



y is a real matrix of transition probabilities. It satisfies Vat = 0. The Hamiltonian for the quantum system is H: 
it is a function of the sequence of states of the classical environment, and is therefore time-dependent. For a fixed 
sequence S the density matrix evolves according to the Von Neumann equation 



^ = -^[HiS),ps] 



(3) 



in units with h = 1. However, we are interested in the density matrix averaged over all sequences. We shall denote 
averages over S by an overbar, so the actual density matrix is p ~ pj- We shall treat both the quantum system and 
the classical environment as finite-dimensional. 

Since a is a classical random variable, this is a classical noise model. The model applies when the noise sources are 
more strongly coupled to an external bath than to the qubit, so that back action of the qubit on the noise sources is 
negligible. The Hamiltonian _ff is a function of a, the state of the classical system, but Vab is independent of p. This 
implies that quantum information that leaves the qubit leaves forever. The conditions under which such a model is 
appropriate have been considered in more detail by Galperin et al. 19]. 



A. Transfer Matrix for a Fixed Noise Sequence 



We wish to compute the qubit density matrix ps (t), given ps (0), for a fixed sequence S. Our first step is to rewrite 
this in terms of the evolution of a generalized Bloch vector rii (t) : 



Ps{t) 



N„ 



1+ J2 n, {S, t) A, 



i=l 



(4) 



where rii is a set of — 1 real numbers, / is the Ng x Nq unit matrix and Xi are the generators of SU (N). The A^ 
are time-independent Nq x Nq matrices and they are chosen to satisfy 



Tr Xi =0, X] ^ Xi, and — Tr AjA^ = 6i. 



(5) 



The Xi form an orthonormal basis for the quantum state space of density matrices under the inner product {pi, P2) = 
^Tr (pip2)- The fact that the A^ are traceless, together with Eq. [4], immediately implies the conservation of probability: 
Trp = 1. 

Consider a short time interval At in which H is constant and the environment is in a fixed state a. The formal 
solution to Eq. [3] is 



Ps (a. At) = U {a, At) p (0) U'' (a. At) . 
with U (a, t) — exp [—itH (a)]. In terms of the A;, Eq. [6] is 

^lafi + ^rii {a, At) Aj,a/3 = Uaj (a. At) C/]^ (a. At) ^I^s + i^rii (0) Xi,js 



(6) 
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where we have temporarily included Greek subscripts for clarity. These indices denote components in the Hilbert 
space of the quantum system. They take on the values a = 1,2,..., A^^. The Roman subscripts take on the values 
i = 1,2, — 1. Both are subject to a summation convention. The identity matrix term cancels out {lap has no 
dynamics) and we have 

Tii (a, At) \^af3 Ua^ {a, At) [a. At) Ui (0) A^,^^. 

We may extract the components of n by multiplying this equation by the i-th generator and taking the trace over the 
Greek indices. Using the trace identity from Eq. Owe find 

Hi (a. At) = ^Uaj (a, At) (a. At) Uj (0) Xj^^s \i3a- 

This is conveniently written as 

(a. At) = T,j {a. At) (0) , 

where 

r„- (a. At) = ^Tr [\,U {a. At) \,U^ {a. At)] . (7) 

Tij (a. At) is the quantum dynamical map (sometimes referred to as the Liouvillian) for the interval At in an envi- 
ronment in state a (Some properties of T are given in App. |A|. From now on the Greek indices will be suppressed; 
operations in the quantum Hilbert space of operators are indicated by matrix multiplication and the trace. 
Thus for the whole sequence S, 

ns{t) = T(ajv) • • • T(a2) T{ai) n(0), (8) 
where t = NAt and Ofe labels the state of the classical environment in the time interval kAt. 



B. Averaging over All Noise Sequences 



We are interested in the generalized Bloch vector averaged over all possible sequences, i.e. n(t) = ns{t) = T{t)n{0). 

To compute T(t), we note each noise sequence S = {ai, a2, . . . , otv} is associated with probability 
Wao,aiVVai,a2 ' ' ' W^ajv_i.aNJ whcrc Wa^.ak+i IS the transition probability of the classical system going from state Uk to 
ak+i at the end of the fc'th time interval, see Fig. [l]and the infinitesimal expansion of W{At) is given by 

W{At) =Ic + VAt (9) 

where Ic is the Nc x Nc unit matrix. Wag^ai is put in by hand for later on convenience. As — > oo it will not 
introduce any errors. 

Let us defined a [NdN^ — 1)] x [Nc{N^ — 1)] tensor T whose element is 

T{ar, ar-l) = Wa^^a,.^, (g> T{ar). (10) 

The averaged Bloch vector is thus given by 

n[t)= J2 ^ai(0)r(a^,ajv-i)---r(ai,ao)n(0) (11) 

ai,a2,...,aiv 

Equivalently, we can utilize the tensor nature of T and write 

n(t) = (x/|r^lv)n(0) (12) 

where (a;/| = [1, 1, . . . , 1] and — [pi{0),p2{0), ■ ■ ■ ,Pn^{0)] act on the clssical environment. This contraction 
amounts to averaging over all the Nc x Nc blocks of , each of which corresponds to the family of time evolutions 
caused by noise sequences that start from ai and end with qn- In this formalism, it is possible to put the classical 
environment into an arbitrarily initial state \if). However, in almost all cases of physical interest, both the initial and 
final states of the environment will be the stationary distribution \ps)- the right eigenvector of V corresponding to 
the eigenvalue zero [l7|. 
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FIG. 1: A fixed sequence S — {2, 3, 1, 1, m, ... , k, 4} of the classical noise. Pi{0) are the initial probability distributions of the 
Nc states. Wai,ai^i are the infinitesimal transition probability between state at and a^+i at the instant iAt. 



This repeated matrix multiplication structure is the key to the formalism. Expand the matrix F as 

F(A<) w / - iHqAt (13) 
so that in the limit Ai ^ with t = NAt held fixed we have 

F^(t) = (/ - iH^At)^ = cxp {-iHqt) , (14) 

and 



n{t)^T{t)n{Q)^{xf\e-^"^'\tf)n{Q). (15) 
Hq is the time-independent " quasi-Hamiltonian" given by 

H„ = I hm ^^^^J ~ ^ = z^T{t = 0). (16) 
^ At^o At dt ^ ^ ' 

It completely characterizes the evolution of the open quantum system. Hq is pure imaginary and non-Hermitian. 
The idea of a non-Hermitian Hamiltonian to characterize dissipation or absorption in quantum systems is not new, 
going back at least to the optical model of the nucleus [28]. However, the implementation has historically been 
phenomenological; the results here are exact. 
We write the eigendecomposition of Hq as 

Hq^Y. ^V- (V'l , 

with \ip) and (-01 being the right and left eigenvectors of Hq. Note that since Hq is not Hermitian, and (-01 are not 
dual to each other as in ordinary quantum mechanics. 

Note lip) is a state of the combined environment-qubit system and the total evolution is given by 

e-''H->' =J2We-''^^' {^\, (17) 
i> 

so that Re((jj^) gives the oscillation frequencies and — Im(a;^) gives the decay rates of the combined system. ip = 
1, 2, [Nc [Nq — l)] . Included in this list of — Im(a;^)'s are the rates for the environment. 
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This formalism provides a means of calculating the dissipative evolution of any quantum system evolving in the 
presence of a classical stochastic process. Furthermore it is completely linear, which means that all of the powerful 
techniques of linear algebra can be brought to bear, including well-controlled perturbation theory. In favorable cases 
such as the one to be considered next, the quasi-Hamiltonian is similar to Hamiltonians familiar from other problems 
in classical or quantum mechanics. The arsenal of methods developed for these situations can then be brought to 
bear. 

We note that the derivation of the quasi-Hamiltonian is similar to the time-slice derivation of the path integral 
approach to open quantum systems, which leads to the Feynman- Vernon formulas for the influence functional 29]. 
The difference in starting points is that the environment here is taken as classical. The difference in end results 
is quite startling, since a non-Hermitian Hamiltonian does not seem to emerge naturally from the Feynman- Vernon 
approach in any obvious limit. 

We note finally that the present formalism is also ideal for the investigation of quantum control schemes such as 
pulsing a qubit. We only need to sandwich the pulsing operators between the evolution operators. Let the pulsing 
operator be a unitary operator R that acts at the time tp with < ip < t. Its action on the generalized Bloch vector 
is given by the [N"^ ^ l) ^ {^q ^ l) raatrix 

R,.^. = iTr X.RXjR^ (18) 

and writing Up — Ic ^R, we have 

n (t) = {xf\ e-»^''(*-*p)Upe-*^«*'' \if) n (0) . 
The generalization to more complicated pulsing schemes is immediate. 



III. QUBIT DEPHASING BY TWO LEVEL SYSTEMS 



A. Quasi-Hamiltonian 



We now proceed to solve exactly the problem of the evolution of the density matrix of a qubit in the presence 
of an environment of AI independently fluctuating TLS that dephase the qubit. Even this simple case of Nq = 2 
and Nc — 2^^ is of great experimental interest. From the exact formulas we will derive qualitative information by 
extracting and analyzing asymptotic expressions in various limits. 

For M statistically independent TLS we can describe the state of the environment by variables s„ (t) = ±1 that 
switch at random intervals at an average rate 7„. n — 1,2,..., A/. The most general expression for the flipping 
probability matrix is W = W\ ® ■ ■ ■ ® Wm with 

or, in index notation (Wn)_^_i = 1 — Pn — ^n, etc. This states that the probability of starting and finishing the 
interval in the -1-1 state is 1 — Pn — <5ri, the probability of starting in the -1-1 state and ending in the —1 state is p + 
etc. We can then write Wn = 1 — Pn + PnTnx — — iSrny, where r„i are the Pauli matrices that act in the state 
space of fluctuator n. The switching rate is 7„ = p„/At, while 6n controls the average occupation of the states. We 
shall focus on the case Sn = 0, the unbiased fluctuators, when we have 

Wn^il-Pn)I+PnTnx- (20) 

and the stationary state is \ps) = [1,1,..., l]/2*^, which is the unbiased fluctuator. 

Here Nq = 2 and the generators of SU{2) are the Pauli matrices ax,(Jy, and CTz. The Hamiltonian of the quantum 
system is 

M 

I I 

H\ 



(t) = ~^Bq(T^ ^ ^ X! 3n <Jz- (21) 

71=1 

This is the case of pure dephasing noise. The more general case H = —Boaz/2 — J2n {t) 9 ■ i7/2 can also be treated 
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by the same method [l^ . Using Eq. [71 we have the 3x3 matrix 



T,,({s4,A0 = -Tr 



CTie'^-^''+^ s„g„)f .At/2^^.g-i(Bo+E s„g„)cr. At/2 



exp 



iLz I ^0 + ^ Sn9n 1 Ai 



(22) 



where is the usual angular momentum matrix: {Lz)^y — —i, iLz)yx = i, and all other = 0. Eq. [22] can 

be derived by direct calculation or by noting that eyi]3 (ia z9 / 2) is a rotation by 9 about the z-axis in spin space. 
Substituting Eqs. and into Eqs. [TUland then using Eq. [THIwe have the quasi-Hamiltonian for this problem: 



M 



ffq = ^ [-iln + ilnTnx - gnTnzLz) - B^Lz 



(23) 



Note s„ in Eq. [22|are replaced by t„z due to the first order expansion. 



B. Single Fluctuator 

We first consider the case of a single TLS, so that M — 1 and Nc — 2. This simple case illustrates all the essential 
mathematical features of the method and the generalization to many independent TLS is almost immediate. We now 
have 

Hq = + i^T^ - gTzLz - B^Lz- 

The problem of qubit evolution has been reduced to the diagonalization of the 6x6 matrix Hq. This is much simplified 
by the fact that [Lz, Hq] = so the problem reduces to a set of 3 smaller problems for Lz = 0, ±1. In these smaller 
problems no manipulations more complicated than diagonalizing a 2 x 2 matrix are required. We treat the smaller 
blocks in turn. 

1. Lz = 0- The 2x2 block of the quasi-Hamiltonian Hq is 

Hq [Lz = 0) = (-17 + ijT^) 
There are 2 eigenvalues and right eigenfunctions that satisfy 

Hq {Lz = 0) l^L^^n) - ^L.=o |«'L.=o) • 
We label them by s = ±1. The right eigenfunctions and eigenvalues are 

(?) 

iUL,=o (s) = -«7 + is"f 

2. Lz = +1. The 2x2 block is 

Hq {Lz = 1) = -i-f + hT.^ - 9Tz - Bq (24) 
The right eigenfunctions and eigenvalues for s — ±1 are 

u;l,=i {s) = -47 - Bo + s\/ - 7^. 

Here C^^ = 2s\J — 7^ (g + s\J — 7^). The corresponding left eigenvectors (5'l,=i| is simply the transpose of 
|*L,=i(s)) and (*L,=i (s) (s')) = -5,,,.. 
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3. 



Hq {L^ = 1) = -il + i"fT^. + gTz + Bq 

Comparison to Eq. [Ml shows that the eigenvalues and eigenvectors for = — 1 are obtained from the Lz — 1 case 
by the substitutions (l/V2, i/V2, 0) (1/^2,-1/^2,0), g„ -gn, and Bq -Bq. 

We now perform the average over initial states and sum over final states. We assume that Pg {t — 0) has its steady 
state values Pi = P_i = 1/2, or \if) = [1, l]/2. However, the average and sum are more conveniently performed by 
taking a partial inner product with the state | of the classical system 



(25) 



which projects onto the quantum subspace. The final 3x3 evolution matrix in this subspace is: 

Note this equivalence between (y^j • |-\/p^) (-^/l ' I*/) '^^^ only be established if the initial states are uniformly 
distributed. 

Using the eigendecomposition of Hq we have 



m = ^ (^1 (,s)) e 



(26) 



so we need to do a sum of 6 terms for each member of the matrix. We now compute each member of {t) in turn. 
For Tzz(t) only the = 0, s = 1 term contributes and we find 

Tzz{t) = cxp [-iujL,=a (s = 1) t] 
^ 1, 

This is simply the obvious statement that for pure dephasing noise, does not decay. In terms of the standard 
relaxation times, this says that Ti = oo. 
T,z{t) = Tz4t)= Tyz{t) = Tzyit) ^ 0. 

For Txxit) and Tyy{t) only the Lz — ±1 blocks contribute. After a straightforward calculation one finds 



Txx{t) = Tyy{t) = exp (-7i) cos B^t 

and finally 

Txy{t) = -Tyxit) = exp (-7t) sinBoi 
So, for example, if n (t = 0) = (1, 0, 0), then we have 
n {t) = (cos Bot, sin Bot, 0) exp (— 7^) 



cos V5 ^ 7 H + 



7 



sin ( — 7^i 



cos [y g - -1 t] + 



cos (yp~^7^i) + 



vV-7^ 



sin ( yff^-^t 



VP~7^ 



sm I V 5 7^ 



(27) 



The first factor is the uniform precession, and the rest of the expression gives the decay and non-uniform precession 
due to the TLS. 

In the analysis below, it will be convenient to deal with the relaxation function F {t) defined by 



so 



F {t) = 7t - In 



Txx{t)^ cos Bot exp[-F(t)], 



cos (Vg^-l^t) + I J , sin ( V.g^ - 7^0 

^ ' yg -Y^ ' 



(28) 



(29) 
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1. Weak Coupling 

This is the case 7 > Note that weak coupUng (smaU g) is the same thing as fast switching (large 7). The 
arguments of the trigonometric functions are imaginary and T{t) is better written in terms of hyperbohc functions: 



Txx{t) = exp (-7i) cos Bot 



cosh {^■\fi^~~<ftj + 



sinh ( \/ — g^t 



and the behavior at long times » 1 /j) is given by 



exp 



(-7 + t 



(30) 



Thus the dephasing rate is 



1 



7- -g^ 



For the extreme weak coupling case 7 ^ 5 we find 



(31) 



1 



5_ 

27' 



which is the standard result from perturbation (Redfield) theory. 
In the short time limit f <C I/7 we have 



(32) 



(33) 



This reuslt is interesting: it shows that the envelope function initially decays quadratically even for a single fluctuator. 
We shall see below that this behavior is completely generic. 



2. Strong Coupling 



When 7 < 5, one has 

Txx = exp (-7*) cosBqI 



cos (^\f(f-~^tj + 



sin I \/ g^ — 7^t 



and at short times t <C 1/ \/ g^ ~ 7^ we have 

Txx = cos BqI 

and 



T{t) = \gH^~\^gh' + 0{t^) 



(34) 



(35) 



(36) 



Note that when the coupling constant g is increased past 7 the relaxation rate saturates at 7. At the same point the 
oscillation frequency bifurcates into the two frequencies Bq ± gj \J g'^ — 7^. 
We stress that Eq. [34] gives a result that is exact at strong coupling. 



C. Many Fluctuators 



The quasi-Hamiltonian is given by Eq. 1231 Again we have [Lz,-ffq] — 0. The quasi-Hamiltonian for each value of 
hz is a sum of operators acting on the individual fluctuators, which is a sign of the fact that they are statistically 
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independent: they do not "interact" with one another. Thus the generahzation from the single fluctuator case is 
almost immediate. We have 



= 1, 



which is a sign of pure dephasing, and 

/ M 



M 



Tcoxit) = Tyy{t) = exp - ^ 7„< cosBoi H 

\ n=l / n=l 

CM \ M 

~-^lnt\ sin Sot Yi 
n=l J n=l 



[Vgl-lnt) + I 2 sin (\/5^-7n*) 

V 9n In 



COS V5„ -l„t) + 



sm V5n - 7„i 



(37) 
(38) 



J. VKeafc Coupling 



If < 7„ for all n, then 



M 



M 



r(0 = E^"*-Ei'^ 



n=l n— 1 

If t ^ 1/ min„ (7„), then we have 



cosh (y-^ll-gli^ 



■ sinh (ys/ll-glt^ 



M 

r(t)^^(7„-v/^F^)t, 



n=l 



SO that the decay is exponential at long times. For extreme weak coupling Qn -C 7n for all n then the Redfield result 
holds: 



1 

^ = i V ^ 
72 2 ^ 7„ ■ 

n=l 



At short times t <C 1/ max„ (7„) 



M ^ M 
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We get deviations from the quadratic behavior at times of order 



t ■ 



EM 2 ' 

71=1 InQn 



(39) 



(40) 



(41) 



Thus the dephasing behavior is essentially exponential rather than quadratic in the weak coupling region. 



2. Strong Coupling 



If 5n > In for all n, then we find 



M 



M 



m = Y.^nt-Y.\n 



n=l n=l 



cos (ys/gl-ilt^ 



\/9l-ll 



sin (ys/gl - lit) 



(42) 



This equation is exact and repre sents a n ew result for many strong-coupling fluctuators. 
At short times t ^1/ max„ sj g^ — l1 



^ M ^ M 



(43) 



n=l 



n=l 
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We get deviations from the initial quadratic behavior at times of order 



2^n=l Indn 



When t > ts the behavior is more comphcated. We write 

M M 



n=l n—1 

where A„ = \/ " In ^^^'^ need to evaluate the expression 

AI M M 



2 2z .Xr) 



^ ^ ri=l Ti=l 



(44) 



(45) 



(46) 



1/2 

with r„ = 1 - i'jn/Ki = kn| exp (iein) ; |r„| = (l + Jn/Xl) and 6'„ = tan^^ {-In/K) 



To this end we note that X]n=i ■^n (-^n^ + ^n) is the result of a random walk with a large number of stpes M. In 



the long-time limit < 3> = 1/ min„ g„ the central limit theorem gives 



M 



M 



n=l 



exp ■ 



n=l 



E 



5^ 



(47) 



(48) 



This expresssion of course neglects any Poincare recurrences that can occur whenever M is finite. In terms of F we 
have 



M 



n—1 n—1 

SO that there is a regime of Gaussian decay at long times when the coupling is strong. 



M 



(49) 



D. Broad-spectrum Noise 



Finally we consider the case when the noise does not satisfy either of the inequalities g„ ^ 7„ for all n. The 
fluctuators can still be divided into Mf fast (weak coupling, g„ < 7„) fluctuators and Mg slow (strong coupling, 
g-m > 7m) fluctuators. As we have seen, this is not a qualitative distinction - rather it corresponds to a change in the 
analytic behavior of the eigenvalues. This does not spoil the solvability of the model. We have 



Mf 



Mf 



n=l 



cosh V7,^i - git 



7n 



M, 



M, 



rn — 1 



m—1 



sinh [^/ll- git 



[Vgl^^^t] + 



Vg-. 



2 _ ^2 



and in the short time limit t <^ts 



Mf 



Ms 



while in the long time limit t^tf. 



r (0 « i E - /^I^) + y E 3^ 



sin (Vgli- lilt 



(50) 



(51) 



(52) 



and the slow fluctuators will dominate at long times. This result is consistent with those obtained in Refs. [20| and 
[l9| in more specific models. 
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IV. COMPARISON TO APPROXIMATE SOLUTIONS 



The most usual way to characterize noise is by its power spectrum. Since our noise sources satisfy Poisson statistics 
and are independent of each other, we have 



Sm (t) Sn (f) = 5jnn exp (-27„ \t - t' \) , (53) 

and the time auto-correlation function of the noise is 

Y^bmz {t)bn, it') = 9^sm {t) s„ {t') = ff' ^ exp (-27„ \t -t'\). (54) 

mn mn n 

The power spectrum is obtained by taking the Fourier transform: 

1 



S{lo)^— dt (t) b, (0) e-^'"*. (55) 

ZTT 



For our case this is 

n=l 

each individual fluctuator follows Poisson statistics and has a Lorentzian power spectrum. Sd (w) is an even function 
of frequency; this is probably the main limitation of our classical model, as quantum noise is asymmetric in frequency 
at low temperatures (26| . In the continuum limit, we find 

Sci{uj)^- dg / d-ip{g,^)—— — (57) 
71" Jo Jo + ^ 

where p{g,j) is the distribution of couplings and rates, defined as 

M 

P{9,l) = YH9 - 9n) In) ■ (58) 

n=l 

When many fluctuators are superposed, we can obtain an arbitrary power spectrum by the proper choice of p (5,7). 
Indeed, even choosing gn — go independent of n we have 

1 2^ 
Sd{.) = -glj^ d7W7)^:^ (59) 

1 



TT 



gl / dr cos cot / dj p {-/) e'^'^K (60) 



Defining p (7) by p (g, 7) — S {g — go) p (7), this equation shows that to obtain p (7) given Sd (i-o), we first invert a 
Fourier cosine transform to obtain the original time auto-correlation function and then p (7) is proportional to the 
the inverse Laplace transform of that. We conclude that as long as the only characterization of the noise is its power 
spectrum, then the results given above provide an exact solution for any S (lu). 

We have already commented on the relation of the present solution to perturbation (Redfield) theory. The exact 
solution agrees with the perturbative results when gn <C 7n for all n and t 3> 1/ min„7„. 

A more interesting approximation is the Gaussian approximation: 



TG{t) = -l Sd (c.) °"\r/; da.. (61) 



r o f . sin^/2 

— / Sd (to) 

2 J -00 {Lot/2) 

See, e.g., Ref. for a derivation. This approximation is valid when noise cumulants of third and higher order vanish. 
For RTNs, this is not the case - there are cumulants of all orders. For a calculation of some of these cunmulants, 
see Ref. [27|. Cumulants of order n for a single noise source are proportional to g", so we expect that the Gaussian 
approximation will break down for large g. Qualitatively, the behavior of Tq (t) may be obtained by observing that 
the function sin^ {ujt/2) / {ujt/2)^ acts largely as a filter function that passes frequencies w < 1/i, so 

/•I/* 

Fg (t) « / Sd (w) duj. (62) 
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FIG. 2: (Color online) Envelope of the free induction decay signal exp [— F (t)] for weak coupling: g=0.01,7 = 0.2. The exact 
result is given by the blue solid line, Gaussian approximation by the green dashed line and the Redfield approximation by the 
red dashed line. The inset gives the short-time behavior. 



Furthermore, the total noise power is proportional to Sd (uj) dio. For this integral to converge (pathological cases 
apart) there must be an upper {ujuv) cutoff frequency for Sd (uj) and a lower (ujir) frequency at which 5c; {lu) rolls 
over and becomes a constant Sd (0). Hence the asymptotic behaviors of Tq (t) are 

There is an initial quadratic decrease of the signal and pure exponential behavior at very long times. 
It is now of interest to compare Tq (t) with the exact F (t) for some interesting distributions p{g,"f)- 
For a single fluctuator, we have 

Sd{cj)^-9 \ , (64) 



so that the Gaussian approximation is [l£ 



1 



[t) = / , . .2 , , ,„,2 (65) 



474 



+ 27< - 1) . (66) 



Our results may now be compared with Redfield theory and the Gaussian approximation (Eqs. [29| . We give the decay 
function exp [— F [t)] in Fig. [2] and [31 For weak coupling g < 7 [Fig. [2] Redfield theory works except at short times 
t <^ 1/7, while the Gaussian approximation is excellent at all times. For strong coupling 5 > 7 the situation is more 
complicated [Fig. [3]. The exact solution develops oscillations that are not present in the approximate solution. Again, 
Redfield theory is poor at short times, while the Gaussian approximation is very good at these times, as already noted 
by other authors [24|. At longer times both approximate solutions have little resemblance to the exact solution. We 
summarize the situation in Fig. 3] The areas of agreement (to within 1%) are given by the white regions. Notice 
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FIG. 3: (Color online) Envelope of the free induction decay signal exp [— F (t)] for strong coupling: g=0.01,7 — 0.002. The 
exact result is given by the blue solid line, Gaussian approximation by the green dashed line and the Redfield approximation 
by the red dashed line. The inset gives the short-time behavior. 



that the normahzation is relative to the initial value of the signal. At long times, the ratio of the exact results and 
the Gaussian approximation can be much different from unity; however, the absolute value of the signal is small and 
may be difficult to observe. It is interesting to note that the discrepancy between the approximate theory and exact 
theory is oscillatory and is not well characterized as a " plateau" . This phenomenon is more closely analyzed in Ref. 

V. CONCLUSIONS 

We have given a general formalism for the dissipative dynamics of an arbitrary quantum system in the presence 
of a classical stochastic process. It is applicable to a very wide range of physical systems. This method has several 
virtues. It is linear, and the close analogy to Hamiltonian systems opens up a toolbox of well-developed methods such 
as perturbation theory and mean-field theory. We applied the method to the problem of a single qubit in the presence 
of TLS that give rise to pure dephasing 1/f noise and solved this problem exactly. This has been done before by the 
method of stochastic differential equations [T9| . However, that method depends on a non- linear parameterization of 
the density matrix that is difficult to generalize. We anticipate that the method can be applied to other quantum 
systems, such as an array of qubits, and also other kinds of noise. 
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FIG. 4: (Color online) Region of validity for Refield and Gaussian approximations. Envelope function exp[— r(f)] calculated 
from Gaussian and Redfield approximates are compared with the exact envelope. For this figure, we set 7 = 1. In the green 
and yellow region, the error of Gaussian approximation is smaller than 0.01. In the red and yellow region, the error of Redfield 
approximation is smaller than 0.01. The blue solid line is where the exact envelope dies off (smaller than 0.01). The horizontal 
black dashed line separates the strong (top) and weak (bottom) coupling region. The red spike in the middle is due to the 
fact that the envelope calculated from Gaussian approximation crosses the exact envelope. The yellow island on top is due 
to the fact that both approximations die off much faster than the exact solution and the exact solution crosses zero multiple 
times, as seen in FigO Note in FigO 3/7 = 5, there will be 4 islands for Gaussian approximation and 5 islands for Redfield 
approximation. The blue boundary has zigzag shape in the strong coupling region. This is due to the oscillations of the exact 
envelope in the strong coupling region. 

APPENDIX A: PROPERTIES OF T AND T 

In this appendix we derive two properties of Tij . 

1. T is a real matrix. This is shown as follows. 

7^* =iTrfT*i7V;(C/t)* 

= iTr.f(C/-)^.JC/- 

= ^Tr Un,U-^ a, 

= ^Tr n,Ua,U^ 
— T ■ 

2. T is an orthogonal matrix. This is proved most simply by noting that that the set of 2 x 2 Hermitian traceless 
matrices Ai form a 3-dimensional real Hilbert space with inner product [Ai^Aj) ~ ^Tr{AiAj). The ai are a complete 
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orthonormal basis for this space. Here is the proof. 

TijTkj = ^Tr {(j.UajU^) Tr {akUajU^) 

= {U^a.U.Uj) {U^akU,aj) 

= ^TV (U^aiUU^akU) 

= ^Tr ((7j(Tfc) 
= Sik- 

The averaged quantity Tij (t) is real since the averaging is over real weights, but Tij {t) is orthogonal only in trivial 
cases. 
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